# Load Necessary Packages
library(foreign)
library(ggplot2)
library(lme4)
library(marginaleffects)

# Load Dataset
# setwd("")
besip <- read.csv("BESIP Long.dta",header=T,sep=",")

# Run Model
besip$monarchwave <- ifelse(besip$monarchwave=="monarchism25",1,0)
besip$cohort_m <- besip$cohort-1926
m.besip <- lmer(monarchism~(cohort_m*monarchwave)+(1|id),data=besip)

# Create Figure 5
## Make Predictions
besip.predictions <- predictions(m.besip,newdata=datagrid(cohort_m=unique,monarchwave=unique))
besip.predictions$cohort_m <- besip.predictions$cohort_m+1926
besip.predictions$monarchwave <- factor(besip.predictions$monarchwave)
## Create Plot
cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
f5 <- ggplot(besip.predictions,aes(x=cohort_m,y=estimate,group=monarchwave,color=monarchwave)) + geom_line(lwd=1.5) + geom_ribbon(aes(ymin=conf.low,ymax=conf.high),alpha=0.2) + theme_classic() + scale_color_manual(name="Year",values=c("#E69F00", "#56B4E9"),labels=c("2016","2023")) + ylim(0,1) + xlab("Cohort") + ylab("Support for Monarchy")
# ggsave("Figure 5.png",f5,device="png",width=6,height=4,units="in")
